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We derive an analytical description of neutrino oscillations in matter based on the Magnus ex- 
ponential representation of the time evolution operator. Our approach is valid in a wide range of 
the neutrino energies and properly accounts for the modifications that the respective probability 
transitions sufi^er when neutrinos originated in difi'erent sources traverse the Earth. The present 
approximation considerably improves over other perturbative treatments existing in the current lit- 
erature. Furthermore, the analytical expressions derived inside the Magnus framework are remark- 
ably simple, which facilitates their practical use. When applied to the calculation of the day-night 
QQ ■ asymmetry in the solar neutrino flux our result reproduces the numerical calculation with an accu- 

' racy better than 1% for the first order approximation. When the approximation is extended to the 

, second order, the accuracy of the method is further improved by almost one order of magnitude, 

^S) ■ and it is still better than 5% even for neutrino energies as large as 100 MeV. In the GeV regime 

characteristic of atmospheric and accelerator neutrinos this accuracy is complemented by a good 
O ' reproduction of the position of the maxima in the fiavor transition probabilities. 
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I. INTRODUCTION 



Neutrino physics has experienced a spectacular progress in the last decade. Many experiments with neutrinos from 
different natural and artificial sources have provided convincing evidence on the existence of neutrino oscillations, 
^ [ a remarkable quantum interference phenomenon taking place at macroscopic distance. Experimental results can be 
'— 'i satisfactorily accommodated within a scheme where at least two neutrinos are massive and there exist a leptonic 
' mixing analogous to the one in the quark sector. From the present data set, two neutrino mass-squared differences 
, and two mixing angles have been determined [l|: {Sm^i = — ml w 8.0 x 10^^eV^,6'i2 ~ 35°) driving solar and 
ly-j reactor neutrino oscillations and (|(5m|2| ~ 2.5 x 10~^eV^,6'23 ~ 45°) which drives atmospheric and long baseline 
, neutrino oscillations. The third angle ^13 and the CP-violating phase remain undetermined. The determination of 
T-H ■ these parameters, as well as the determination of the neutrino mass hierarchy, will be the main goals of the next 
\ generation of experiments. We are thus entering into a new stage characterized by high precision measurements. 

■ In turn, the interpretation of the forthcoming results will require more careful theoretical descriptions of neutrino 

■ oscillations that incorporate sub- leading processes. 

00 A subject of particular interest within this context, refers to the matter effects on the flavor transformations for 

neutrinos propagating through the Earth. The problem has been investigated by direct numerical integration of the 
^ , equation that governs flavor evolution in a medium. Yet, analytic calculations have been implemented to simplify the 
numerical computations greatly and also to gain a better understanding of the underlying physics. Many of these 
studies have been carried out under the assumption of one or several layers of constant density. Extensions for a 
varying density have been developed on the basis of the perturbation theory for oscillations, both in the low 0, [1, 0] 
energy and high energy regime [5, 6, ■F'^r low energy neutrinos the perturbative solutions were found in the basis 

of the mass eigenstates, while in the high energy limit the method was formulated in the flavor basis. In this work ,we 
present a novel analytic description of the effect based on the Magnus exponential expansion of the time-displacement 
operator to), which make possible a unifled treatment of the problem and give us precise simple formulas for both 
energy ranges. 

The evolution of the flavor amplitudes of a neutrino system may be conveniently described in terms of the operator 
U, which satisfles the Schrodinger-like equation [§] 

ih^it,to) ^ H{t)Uit,to), (1) 
at 
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with the initial condition U{to, to) = I. Later we shall give the explicit form for the matrix H{t) in the MSW theory. 
The Magnus expansion [1^] supplies a method for finding a true exponential solution of Eq. (1) of the formZ^ = exp(ri) 
(i.e., without time ordering). The operator Q satisfies its own differential equation which in turn is solved through a 
series expansion; fl — ^"j where r2„ is of order h^". The first two terms are explicitly given by 

= - f dt'H{t'), 

^Jto 

= "i/^*'/ dt"m'),Hit")]. (2) 

Due to the antihcrmitian character of every operator r2„, truncating the series for Q at any order gives a unitary 
approximation toU. This is briefly what we shall need to know about the Magnus expansion for its present application; 
further details about the formalism and recursive procedures for building up the successive terms can be found in the 
specific literature [111, IHi ■ 

Here, we use the first and second order Magnus approximation to seek solutions to the problem of 2v oscillations 
in a medium with an arbitrary density profile, which is symmetric with respect to the middle point of the neutrino 
trajectory. The method is based on a formalism that was developed several years ago in order to incorporate non- 
adiabatic effects in the flavor transitions of neutrinos that propagates trough a matter-enhanced oscillation region 
(l^ . The main idea is to follow the time development of the system in the adiabatic basis of the instantaneous energy 
eigenstates and to incorporate the corrections to adiabaticity trough the Magnus expansion. In the Magnus 
approximation was used to deal with the same problem but in the base of the (non-evolving) mass eigenstate. When 
applied to the calculation of the day-night asymmetry for solar neutrinos the method renders a simple formula for the 
regeneration factor, which has a better agreement with numerical calculations than those derived by usi ng p erturbation 
theory. The approach we are presenting now is not only more accurate than the one developed in [l5| . but is also 
valid in a much wider energy interval, allowing for a unified description of the Earth effect on the oscillations of both 
low and high energy neutrinos. 

The paper is organized as follows. In the next section we describe the basic ingredients of the formalism and 
derive the formula for the flavor transition probability in a medium with varying density. In Sec. Ill we present 
two applications of physical interest. In the first one it is shown how the regeneration phenomenon of solar neutrinos 
traversing the Earth can be conveniently accounted for by our present approach. In the second application, we examine 
the influence of the terrestrial matter on the probabilities for <-> f^^T- transitions. Sec. IV contains the conclusions. 



II. FORMALISM 



Typically, the quantity of interest is the probability P^^ of observing an electron neutrino at a distance L ~ 
tf — to from a source {H = c — 1). If |j^(i/)) represents the neutrino state at time tf, then P^^ = \{i'eW{tf))\'^ = 
\{ve\l^{tf,to)\i'{to)}\'^, where \i^{to)) denotes a certain initial state. We consider oscillations between two neutrino 
flavors, let say and Va- In the relativistic limit and after discarding an overall phase, the Hamiltonian of the system 
in the flavor basis {|i^e), Wa)} can be written as 

Ao f -cos2e sm20\ V{t) f I 0\ 

where is the mixing angle in vacuum and we have defined Aq = 5rn? /2E, with E the neutrino energy and 5w? the 
squared mass difference. The effect of the medium is accounted for by means of V ^ the difference of the potential 
energies 14 and Va- In normal matter, to lowest order in the Fermi constant Gp, we have V{t) = Ve{t) — Va{t) = 
\/2 G pneit) , where rig is the number density of electrons along the neutrino path. 

The evolution operator in the flavor basis can be expressed asU{tf,tQ) = Um{tf) ^"^(i/i ^o) Ujj^{to), in terms of the 
corresponding operator U-^{t, to) in the adiabatic basis of the (instantaneous) eigenstates {\i/im{t)), \i'2m{t))} of H{t). 
Here, 

U^it)=( '^"^^-W '^i-'^rl) (4) 

^' \~Sm0m{t) COsBmit) J ^ ' 

is the orthogonal transformation that, at each time, diagonalizes the matrix in Eq. ([3]). The mixing angle in matter 
0„i(t) is given by 

Ao sin 26* , , 

sin2^„(t) = ^-^, (5) 
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where 

A„(<) = Ao - cos 261)2 +sin2 26* (6) 

stands for the difference between the energy eigenvalues and we have introduced the non-dimensional quantity e{t) = 
V{t)/Ao = 2EV{t)/Sm^. 

If V{t) is symmetric with respect to the middle point of the neutrino trajectory t — {tf + to)/2, then 6rn{tf ) = 
Omito) = C and 

U{tf,to)^ Urn (to) U^{tfM) Ul (to). (7) 

This is the situation for the Earth, in which the angle evaluated at the surface. In what follows, we restrict 

ourselves to such a case and find an analytical expression iorU{tf,to) in terms of U-^ {tf,to) calculated by means of the 
first-order Magnus approximation. We follow the procedure presented in Ref. adapted to the present situation. 
To make the work self contained we repeat here some of the steps presented there. 

The evolution operator in the adiabatic basis is a 2 x 2 matrix that obeys Eq. ([1]), with the Hamiltonian 

H^{t)^HD{t)-tUUmrnit), (8) 

where H]j{t) = —^Am{t)az is a diagonal matrix whose elements are the eigenvalues of Eq. ([3]) and U^{t)lJm{t) — 
i0m{t)cry. Here, dot means differentiation with respect to time and and Cj, are Pauli matrices. 

Neglecting the second term in Eq. ([8]) corresponds solving the problem in the adiabatic approximation. In any 
case, the time dependence generated by i?£)(t) can be integrated exactly by a change of the representation, which is 
readily accomplished by means of the unitary transformation ^^■^(i, to) = 'P{t,to)Up{t,tf)), where 



Vit,to) = exp(^-i^ dt'Eoit')^ 



e5<^*o-.t 



(9) 



with 



In the new picture the evolution operator obeys 

I — r— 

dt 

where 



\t'A^{t'). (10) 



H^m:^, (11) 



Thus, we have removed not only the diagonal part, but the remainder of the Hamiltonian gets a simple structure 
which facilitates the algebraic manipulations that follows. 

In general, it is not possible to solve (jlip exactly and one has to rest on some approximation in order to determine 
U^. We employ here, with this purpose, the Magnus expansion and write = e^. Without loss of generality one 
can take detU^p = 1 and, therefore, to any order the Magnus operator has to be of the form = — i ct. ^, where the 
components of vector a are the Pauli matrices and ^a,, ^y, and are real coefficients whose specific forms depend on 
the order of the approximation used to determine f2 in terms of . Consequently, we have 

Z^^ = cose/+^r!, (13) 

where / is the identity matrix and ^ = ^ _|_ Ymia Eqs. (fT3|) and ([9]) it can be shown that U'^ is of the 
general form 
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with the condition \U^\^ + jZij^P = 1. The evolution operator in the flavor basis has the same matrix structure as it 
is easily checked by substituting into Eq. ([7]). 

Subsequently, we put = fli + Q2 and flnd Qi^2 by means of the formulas given in Eq. ([2]) evaluated with the 
Hamiltonian of Eq. (|12p . Proceeding in this manner, and after some algebraic manipulations, we arrive at: 



cos^-zsin£^ e"^'-'/ z sinf^ 

U''{tf,to)= I ^ J ^ / V\ vA h (15) 



i sin^^ ( cose + i sin^^ ) e""^*"**/ 



where ^ = y ?fi) + ^^2)- Non-adiabatic effects on the evolution of the flavor amplitudes are incorporated through the 

quantities ^(i) and ^(2): which come from the flrst-order and the second-order Magnus approximations, respectively. 
They are given by 

fd) = 2J^"di'^sm#^,„ (16) 

The above expression for e^i) was obtained by taking into account that V{t) — V{2t — t) for a potential that is 
symmetric with respect to the middle point of the neutrino trajectory. In this case, 9m{t) ~ — dmi^t — t) and 
dt'Omit') sin 0(^4/ = dt' Omit') siii4>i~>t', while J^^^ dt'9m{t') cos(t>t^t' = 0. By integrating by parts, Eq. (flB]) 
can be rewritten as 

/•*/ 

= 2em{tf) Sincpt^t^ -2 dt' 0rn{t') Am{t') COS (j)t^t,. (18) 



We see that U'^, as approximated by Eq. p5|) . has the form of the general matrix given in Eq. ([T4|) . This guarantee 

that the unitary condition lA^ ^ = U'^'' is verified to second order. As mentioned in the introduction, this is an 
important quality of the Magnus expansion that remains true at every order. In addition, the off-diagonal elements 
of matrix p5|) are purely imaginary, i.e., — ~^i2: but in general this will not be verified when contributions of 
higher order are included. The same considerations apply to matrix U. 

Suppose that |i^(io)) = o:\ve) + P\va), with a and /3 non-negative (real) numbers satisfying + = 1 then, taking 
into account the relations between the Uw {£, £' = e, a) just indicated, we find 

-P.e = a' + (/3' - a^KlmUeaf + 2a(3 {ImUee) {luiUea), (19) 

with 

ImWee = cos26'^ImWi;j -Hsin26'°„ImWi;^ , 

ImZ^ea = - sin 26*55, ImZ^lj -1- cos 26lJ^ImZii^, (20) 



where, according to Eq. ([T 

ImUf^ = cos C sin ^f^^^ - sin ^ ^ cos (j)t-,tf , 

ImU^2 = sinf^. (21) 

As we see, to this order, only the imaginary parts of the matrix elements of the evolution operator are relevant to 
the calculation of P^^ . The result for the lowest-order Magnus approximation is obtained by putting ^(2) = in the 
previous expressions for the imaginary parts oihtyi andZYj^- 

Formula p^ . with the imaginary parts oiUee and Uea given by Eqs. ([^0]) and (HI]), represents our main result. It 
provides an elegant and systematic description of neutrino oscillations in a medium with a symmetric, but otherwise 
arbitrary, density profile, which is valid for a wide range of energies. In order to illustrate its usefulness, in the next 
section we will apply it to two situations of physical interest where the 2v oscillations are suitable to account for the 
leading process: i) the regeneration effect of solar neutrinos when they goes trough the Earth and ii) the oscillations 
of high-energy neutrinos in the Earth. 
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III. APPLICATIONS 



A. Day-Night Neutrino Asymmetry 

The relevant quantity in connection with the solar neutrinos is the probability for a neutrino born as a in the 
interior of the Sun, to remain as a Ve at the Earth. The oscillation parameters controlling the leading effects are 
9 = 9i2 and 6m? = 5m\2 [3- If the phase information is lost, as will typically happen for neutrinos traveling a long 
distance to the detection point, then according to the LMA-MSW solution the averaged survival probability for the 
electron neutrinos can be written as [lB| 

P(l'e V^) = Sil? e + cos 20 C0S^6l° - cos 26g freg-, (22) 

where 9q denotes the matter mixing angle at the production point in the interior or the Sun. The regeneration factor 
freg — ^2e — siu^ 9 represents the terrestrial matter effects expressed as the difference between the probability for V2 
to become Ve after traversing the Earth P2e = P{v2 ^ t^e) ~ |(^'e|^(i/, ^o)|t^2)P and the same probability in vacuum 

\{Ve\y2)? 

We will determine freg by calculating P2e in terms of Eq. p^ . with |i^(<o)) = |f^2) = sin6'|i^e) + cos9\Vf^). Accord- 
ingly, we get 

P2e = sin^ 9 + cos 29 (ImZ^e^)^ + sin 29 IniUee ImUe^ , (23) 

and from this result 

freg = COS 26*;^ COS 26'°, (Im )2 + sin 26i°„ siu 26l°„ (Im )2 

- sm{29l + 2C)(ImZii1)(ImZilt) . (24) 

Here, 9^ = 9^— 9 is the rotation angle that relates the basis of the mass eigenstates ), W2)} with the adiabatic one, 
evaluated on the surface of the Earth. For a constant potential ^ = and, taking into account that sin 29ra = e sin 29m, 
we recover the exact expression for the regeneration factor in a uniform medium 



freg = Sq siu 29^ siu 



^(i/-to) 



(25) 



On the other hand, for the LMA parameters of the solar neutrinos £ ^ 1 within the Earth. In this limit, 

26'„,.(t) = 29 + sin29e{t) + 0{e^), (26) 

as can be shown by using Eq. ^ and 

d9rr,. sin^ 29m. 



de 2 sin 26* ' 

Substituting Eq. ^ into Eq. dH]) we find ^ = -/ + 0(e), where 



(27) 



/ = sin 29 / dt' V{t') cos <j)i-,t' ■ (28) 



In this way, neglecting quantities of 0(e) and higher everywhere, except in the adiabatic face (f>t^t', we arrive at 

freg — sin 21 siu 29 sin (f)t^t + sir? I cos 29, (29) 

which coincides with the expression for freg that was derived in Ref. [Tsj by applying the Magnus approximation to 
solve the equation for the evolution operator in the basis of the mass eigenstates. As pointed out there, by keeping 
the lowest terms of the expansion in /, Eq. ((29)) reduces to the result obtained by means of the perturbation theory. 

In order to compare our results with those corresponding to the first and second order in the e-perturbative 
expansion, we consider a simplified model for the electron density inside the Earth, the so called mantle-core-mantle 
[17| . In this model the electron density is approximated by a step function and the radius of the core and the thickness 
of the mantle are assumed to be half of the Earth's radius i?® : 

5.95 cm-3, r < i?e/2 
ne{r)=NA{ . (30) 

2.48 cm-^, i?e/2 < r < i?® 
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Following Ref. we introduce the function 

^^^^ = T^Tn k-r^ (^) - f^Tr'^ (^)l ' (31) 

freg{E) L « "J 

where frtg^^^ is given by a certain (approximated) analytical expression, frlg"''^^'^ is obtained from the exact (numerical) 
solution and freg{E) = 1/2 Eq sin^ 6* is the average regeneration factor evaluated at the surface layer. Essentially, 5{E) 
represents the relative error of the approximated expression. 

In Fig. [1] we show 5{E) as a function of the energy for neutrinos that cross the Earth through its center. For 
the "solar" oscillation parameters we take 5m\^ = 8 x 10~^ eV^ and tan^ 9i2 = 0.4. As shown there, 5{E) for the 
different Magnus approximations is always smaller than those corresponding to the perturbative calculations. As 
already pointed out in Ref. jl^] , although the error associated with Eq. (|29p increases with energy it remains smaller 
than ~ 2% for i? < 14 eV. The lowest-order adiabatic result derived by doing ^(2) = in Eq. p4)) works even better, 
reducing the relative error to less than 0.5% within the same energy interval. When the calculations in this basis are 
carried out up to the second order, the accuracy improves notably and the error is reduced by almost an order of 
magnitude as compared to the one for the first order formula and remains less the 5% for energies up to 100 MeV. The 
last interval comprises the energy values that are typical for neutrinos originated in supernovae explosions. It should 
be noticed that our treatment works comparatively well in the whole range of energies, whenever the two neutrino 
approximation remains valid, which requires E <C 5m\^/ (2Fsin(?i3 cos6'i3 sin0i2). In order to illustrate this point, in 
Fig. [5] we plot 5{E) for energies as large as 10 GeV, both for the first and second order calculations corresponding to 
the adiabatic Magnus expansion and the perturbative approach. 

B. High Energy Neutrinos 

In this subsection we apply the present formalism to the oscillations of high-energy (E > 1 GeV) neutrinos that 
go across a material medium with a symmetric density profile. If we assume that ^13 is not very small, then the 
quantity 6m\^/2E can be safely discarded in the equation governing the flavor evolution of a 3i^-system In this 
case, the mixing angle 9i2 does not play any role and the problem reduces to an effective one of two states and 
\va) = sin 023 + cos 023 Wt) , whcrc the matter oscillations are driven by the parameters 5m? = 5m\^ and = ^31. 

We focus hereafter in the transition probabilities P{y^J^ v^) — siv?92zP{va ^e) and P{vr —> Ve) — 
cos^ 6*23 P{va — > Ve). Now, \v{tQ)) — \va) and according to Eq. (jl9p . with a = and /3 = 1 we have 

^ {cos2el^luYU^^-sin2dl^luiU({f, (32) 

where ImU^ and ImU^ are determined from Eq. (PT|) . 

Suppose that V S> Aq; then, e 3> 1 and we can implement a perturbative expansion in 1/e for a varying potential. 
Accordingly, 

20,„ ^ TT - i sin 20 (33) 



and 

^ Ao sin 29 J dt' cos ^f^t, + '^(~)- (34) 
Using the last two equations and keeping at most terms of 0(1) in 1/e (except in the phase 0f^f/), Eq. ([32]) becomes 



sin ^Ao sin 29 dt' cos (j)t^t' 



-I 2 

(35) 



in the first-order Magnus approximation (^(2) =0). It is pertinent to note that the perturbative result presented in 

r 1 2 

Ref. [5] P{va Ve) = Ag sin^20 dt' cos 4>t-^t' follows immediately from Eq. (|35p when the sine function is 
replaced by its linear approximation. 

The expression in Eq. (|35p corresponds to the result derived by working directly in the flavor basis, following an 
approach similar to the one we used in [l^. This requires the factorization of the evolution operator as U{t,to) = 
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FIG. 1: The relative error S (Eq. (|3T|) as a function of the energy for a neutrino crossing the center of tlie Earth. The lower 
panels correspond to the the envelopes of \5\, i.e., to the maximum error to be expected at a given energy. The oscillation 
parameters are 5m2i = 8 x 10"^ eV'^ and tan'^ 9i2 = 0.4 and the density profile has been approximated by the core-mantle-core 
model [l^]. (a) and (b) correspond to first and second order of the perturbative approach respectively, (c) corresponds to Eq. 
(|29[) and (d) and (e) correspond to the first and second order Magnus calculation in the adiabatic basis respectively. 



V'' {t,to)U-p{t,tQ), where V is the same diagonal matrix given in Eq. and the determination oi U-p{t,tQ) in terms 
of the lowest-order Magnus approximation ^//^(t, to) = exp[— i J^^ dt' H-p{t')], with the Hamiltonian 

Hvit) = V{t,h)[H{t)- HD{t)]V\tM) 
^ . / e^^'o-t \ 

= ^^^"2%e-^'o^. )■ ^^^^ 

In the above equation, Hd is again a diagonal matrix whose elements are the eigenvalues of Eq. ([3]) and the second 
line has been obtained by using Am(i) = — Aocos26']. Proceeding in this way, the matrix representation for 

U{tf,to) becomes 

cos^(i)e'"^'-'*/ isin^(i) \ , . 

I sm^(i) cosC(i)e ^' V j 

with ^(1) calculated according to Eq. ((34)) . From the last expression we see that ImUae — sin^(i) and formula |35)) 
follows immediately when this result is substituted into P{va t^e) — (ImZYae)^. 

In Fig. [3] we plot P{i'a Vf.) as a function of for the same model of the Earth's density profile used in the 
previous section. We show the numerical calculation together with the analytical approximations corresponding to 
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FIG. 2: The relative error S (Eq. (|3ip l as a function of the energy (up to 10 GeV) for a neutrino crossing the center of the 
Earth. The oscillation parameters are Jm^i = 8 x 10~^ eV^ and tan'^ 9i2 = 0.4 and the density profile has been approximated 
by the core-mantle-core model. Curves plotted in left panel correspond to first and second order of the perturbative approach 
and the ones of right panel correspond to the first and second order Magnus calculation in the adiabatic basis. Note the different 
y-axis scales used in the graphics. 




E [GeV] E [GeV] 

FIG. 3; P{h'a —* i^e) as a function of the energy for a neutrino crossing the Earth passing by its center (left panel) and for 
a trajectory of Nadir angle Q = 26° (cosG = 0.9) (right panel). The oscillation parameters are Sm^i = 2.5 x 10~^ eV^ and 
6^13 = 10°. (a) corresponds to the perturbative approach (see Ref. [^), (b) to Eq. (|35p . (c) and (d) correspond to the Magnus 
approximation implemented in the adiabatic basis for the first and second order respectively and (e) numerical calculation. 
Our approximation reproduces very well both, the value and position of the maxima of the numerical calculation. 

tlie Magnus expansion and to the perturbation theory at first order in 1/e. From tlie figures, it becomes clear that tlie 
formula derived by means of tlie Magnus expansion implemented in tlie adiabatic basis gives a better approximation 
than tlie perturbative method. Moreover, they never give probabilities higher than one, a pathology presented by the 
perturbative expressions as can be seen in the left panel of Fig. [3] This behavior remains true also for the formula 
given in Eq. (pS)) . but in this case the approximation breaks down numerically for energies E Ri {5 — 10) GeV, that 
corresponds to the resonance condition V ~ Ao, for Sm^j^ = 2.5 x 10-3 eV^ The same limitation applies to the 
perturbative result quoted above. 

If the mixing angle 9i3 were vanishingly small or zero, then the problem is also described in terms of an effective 
two-state system {|i^e), |i/f,)}, with — cos023|f^/j) ~ sin(?23 Ij^t)- In this case, the transition probabilities are 
Pii^fi ~^ Ve) = cos^ 62-i Pivb Ve) and P{i't — t^e) — sva^ 623 P{vb t^e), where P{i'b — > Ve) can be computed by 
the same expression given in Eq. (|32p . but with the oscillation parameters Sm? = Srn^i and 9 = 612. From the 
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FIG. 4: P{vfy —f v^) as a function of the energy for a neutrino crossing the Earth passing by its center (left panel) and 
for a trajectory of Nadir angle O = 26° (cosO = 0.9) {right panel). The osciUation parameters are Sm2i = 8 x 10~^ eV^, 
tan'^ 9i2 = 0.4, and 6^13 = 0. (a) corresponds to the perturbative approach, (b) to Eq. (|35p . (c) and (d) correspond to the Magnus 
approximation implemented in the adiabatic basis for the first and second order respectively and (e) numerical calculation. 



curves plotted in Fig. [U it is again evident that the analytical expression derived by means of the adiabatic Magnus 
expansion gives the best approximation to the exact (numerical) result. 



IV. CONCLUSIONS 



We have shown that the Magnus expansion for the evolution operator implemented in the basis of the instantaneous 
energy eigenvalues, provides an elegant and, at the same time, efficient formalism to describe neutrino oscillations in 
a medium with an arbitrarily varying density profile. This approach incorporates in a simple way the Earth matter 
effects on the transition probabilities for neutrinos with a wide interval of energies, making possible a systematic 
description of such effects in the case of solar and atmospheric neutrinos. In both cases, the results are considerably 
more accurate than those derived by different perturbative calculations in the low and high energy regimes. The 
same formalism can be applied without additional difficulties to the study of other situations of physical interest, like 
supernova neutrinos or long baseline experiments with accelerator neutrinos. 
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